Функции построения ЛАЧХ и переходной функции

Clear["Global`*"]

Bode[Φ_,f0_,fn_]:=Module[{ω},
GraphicsColumn[
{
LogLinearPlot[20*Log10[Abs[(Φ/p)/.{p->I*ω}]],{ω,10^f0,10^fn},AxesOrigin->{10^f0,0}],
LogLinearPlot[180/Pi*Arg[(Φ/p)/.{p->I*ω}],{ω,10^f0,10^fn},AxesOrigin->{10^f0,0}],}
]]

Step[Φ_,tn_]:=Module[{t},Plot[InverseLaplaceTransform[Φ/p,p,t],{t,0,tn},PlotRange->Full]]

Продольное движение

(*Вектор X продольного движения*)
Xprod={ ωz[t], α[t], H[t], ϑ[t]}

(*Матрица А продольного движения*)
Aωz={c1+p,-1,0,-1};
Aα={c5 p+c2, p+c4,c6,0};
AH={0,0,p,0};
Aϑ={0,0,-c6,p};
MatrixForm[Aprod=Transpose[{ Aωz, Aα, AH,Aϑ}]]

(*Вектор U продольного движения*)
Uprod={δв[t]}

(*Матрица B продольного движения*)
MatrBprod={-c3,0,0,0,0}

(*Уравнения продольного движения*)
MatrixForm[prod=Table[ Sum[Aprod[[i,j]]*Xprod[[j]],{j,1,4}]==MatrBprod[[i]]*Uprod[[1]],{i,1,4} ]]

Боковое движение

(*Вектор X бокового движения*)
Xbok={ ωx[t],ωy[t],β[t],γ[t],ψ[t]}

(*Матрица А бокового движения*)
Aωx={p+b1,b6,-b7,-1,0};
Aωy={a6,p+a1,-1,0,-1};
Aβ={b2,a2,p+a4,0,0};
Aγ={0,0,-b4,p,0};
Aψ={0,0,0,0,p};
MatrixForm[Abok=Transpose[{Aωx,Aωy,Aβ,Aγ,Aψ}]]

(*Вектор U продольного движения*)
Ubok={δн[t],δэ[t]}

(*Матрица B продольного движения*)
Bδн={-a5,-a3,-a7,0,0};
Bδэ={-b3,-b5,0,0,0};
MatrixForm[Bbok=Transpose[{Bδн,Bδэ}]]

(*Уравнения продольного движения*)
MatrixForm[prod=Table[ 
Sum[Abok[[i,j]]*Xbok[[j]],{j,1,5}]==Sum[Bbok[[i,j]]*Ubok[[j]],{j,1,2}],{i,1,5} ]]

(*Зададимся коэффициентами*)
c1=0.587;
c2=45.8;
c3=13.6;
c4=1.21;
c5=0.204;
c6=5.721;
a1=1.64;
a2=40;
a3=7.56;
a4=0.62;
a5=7.36;
a6=0.985;
a7=0.0278;
(*a8=1;*)
b1=7.235;
b2=133;
b3=22.9;
b4=0.0191;
b5=-1.05;
b6=0.131;
b7=0.004;

(*Оценим возможность изолированного исследования движения "рыскание-скольжение" и движения по крену*)
If[ b1 (a1 a4+a2)≥9 (a1 b2 b7-b2 b6-a6 (a4 b6+a2 b1)),Print["Можно"],Print["Нельзя"]]

(*Коэфициент демпфирования *)
ξβ=(a1+a4)/(2*(a2+a1 a4)^0.5)
If[(ξβ<0.2 || ξβ>0.4),Print["Нужен демпфер рыскания"],Print["Не нужен демпфер рыскания"]]

(*Частота собственных колебаний свободного движения по углу рыскания*)
Tβ=1/(a2+a1 a4)^0.5
Ωβ=1/Tβ

Изолированное движение самолета в углах "рыскание-скольжение"

(*Вектор X бокового движения*)
Xbok1={ ωy[t],β[t]}

(*Матрица А бокового движения*)
Aωy={p+a1,-1};
Aβ={a2,p+a4};
(*Aωy={p+a1,a2};
Aβ={-a8,p+a4};*)
MatrixForm[Abok1=Transpose[{Aωy,Aβ}]]

(*Вектор U продольного движения*)
Ubok1={δн[t]}

(*Матрица B продольного движения*)
Bδн={-a3,-a7};
MatrixForm[Bbok1={Bδн}]

(*Уравнения продольного движения*)
MatrixForm[prod=Table[ Sum[Abok1[[i,j]]*Xbok1[[j]],{j,1,2}]==Bbok1[[1,i]]*Ubok1[[1]],{i,1,2} ]]

(*Запишем передаточную функцию Wy/δн замкнутой системы*)
W0z=Simplify[Det[{Bδн,Aβ}]/Det[{Aωy,Aβ}]]

(*Построим передаточную функцию Wy/δн замкнутой системы*)
Step[W0z,6]

(*Запишем передаточную функцию Wy/δн разомкнутой системы*)
W0r=Simplify[(W0z/(1-W0z))]

(*Построим ЛАЧХ Wy/δн разомкнутой системы*)
Bode[W0r,-1,3]

(*Рассчитаем параметры демпфера рыскания*)
μ=((2* (Sqrt[a2+ a1 a4])-(a1+a4))/a3)

Tn=(a3/(a3 a4-a2 a7))

(*Обратная связь*)
Wu=Simplify[μ*((Tn*p)/(Tn*p+1))]

(*Запишем передаточную скоректированную функцию замкнутой системы Wy/δн*)
W1z=Simplify[(W0z/(1-W0z*Wu))]

(*Построим передаточную замкнутой скоректированной функцию Wy/δн*)
Step[W1z,6]

(*Запишем передаточную функцию Wy/δн разомкнутой скоректированной системы*)
W1r=Simplify[(W1z/(1-W1z))]

(*Построим ЛАЧХ Wy/δн разомкнутой скоректированной системы*)
Bode[W1r,-1,3]